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In the previous semi-annual report, a Fourier analysis method for unsteady 
aerodynamic modeling was presented in the AIAA paper No.91 -2867 (Ref.1 , a copy of the 
version that was accepted by the AIAA Journal being attached.). The results from a set 
of test data for a 70-deg. delta wing reported in Ref.2 as well as two other sets of test 
data for two-dimensional airfoil sections reported in Ref.3 had been used to verify this 
method of modeling responses of harmonic motions at different reduced frequencies. In 
addition, some other cases for a 70-deg. delta wing’s harmonic ramp motions also were 
calculated by indicial formulation in the past work. 

During this reporting period (9/1/91 -2/29/92), two more sets of test data obtained at the 
NASA Langley Research Center have been used to show the generality of the Fourier 
analysis method and validate the indicial formulation of time integration. They are 
discussed below. 

1. Fourier Functional Analysis 

In the present method, the aerodynamic response can be written as 
Cl = Cave 

+ Ei id + E 2 id + Ci * (H-j -j cx + H 2 id ) * (1 - PDi) 

p • • 2 

+ Ei2&2 + ^22^2 + ^2 * (H-|2® + H22®® + ^32® ) 

. (1 - PD 2 ) (1) 

+ Ei 363 + E 23&2 + C 3 * (Hi 3 a 3 + H 2 30 t 2 d + H 33 <xa 2 + H 43 a 3 ) 

* (1 - PD 3 ) 
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From the past experience of analyzing five sets of experimental data, it was found that the 
final values of Ey and Hy could change significantly by using different initial values for Ey 
and Hy as well as different reference values for Cj. This was mainly because that the 
gradient method for finding the best Ey and Hy was unable to locate the global minimum. 
Furthermore, the denominator of the Padd approximant must have negative roots only so 
that in time domain the corresponding exponential terms will die out at large time. This 
requirement can be regarded as a constraint in the optimization process. 

To remedy these two problems and also make the present method more general and 
user-friendly in analyzing any given set of test data, another optimization loop has been 
added outside the existing loops for the purpose of choosing proper initial values for Ey 
and Hy and the best Cj automatically. 

Other task to improve the accuracy of the present method is to include the static test 
data as additional dynamic stall data at a very low reduced frequency, such as k=1 .OE-6. 
The main purpose of this implementation is to avoid possible poor extrapolations at low 
k in the integration of indicial integral. 

Two sets of harmonic motions test data ,one for a 70-deg. delta wing and the other for 
a F-18 model , were analyzed by the updated version of the present method. The results 
are presented in Figs.1 and Figs. 2 respectively. Five Fourier terms are used in this 
analysis. All results were done by the same set of build-in initial data for Ey.Hy and Cj 
without any try-and-error effort by users. The results show that the updated version of the 
present method is able to capture all hysteresis effects. Comparing with the up-stroke 
data , most of the down-stroke data are modeled with less accuracy. The reason is that 
the trend of the hysteresis behavior on down-strokes is not as consistent as those on up- 
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strokes (see Figs.7). This may imply that a higher order Pad§ approximant could be 
needed to model responses which have more complicated hysteresis effects. It is noted 
that the mismatched part of F-18 C m response at k=0.0075 is due to the even more 
inconsistent trend of the hysteresis behavior occurring in the region of high angles of 
attack. 

2. The validity of indicial formulation of time integration 

As indicated in Ref.1, an indicial formulation for any arbitrary motion is written as 


m 

Cl^) - CLjndicJt'-^ <*(*). “Wlx-O + c ave + E ( E 1j«j + E 2j«j) 


j-1 


+ E i^*(1-a 1j e- a 3i( t/ ^-a 2j e- a ^ (t/ ^) SSSL dx 
j.-l^O da 1J 4 dx 

+ J_ E ft' d(a -V l.(l-a 1je -^-aae-^ (t/ ^) ^ dx 
v oo j- 1 *'0 dd 11 dx 


( 2 ) 


To show the validity of the above formulation, three types of motion have been analyzed. 

(1) Harmonic motion and harmonic ramp motion: 

The first type of motion is used to compared with Fourier modeling results; meanwhile, 
the second one is used to show the agreement with the static value at the time when the 
motion stops. 

As it was known in the past study, discontinuity could happen in the calculation of time 
integral if the given motion has a sudden change in a . This can be easily solved by 


slightly increasing the amplitude of original model (a 0 ; see Eq. 3a), such as by 2.5 
degrees. In the following calculations, the amplitude of harmonic model which originally 
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was 30. degree has been set to be 30.5 degree for harmonic ramp motions and 32.5 
degree for constant rate pitching motions respectively. It should be emphasized that this 

does not change the actual instantaneous values of a and a in the indicial time 

integration. It merely changes the values of equivalent frequency (k) and phase angle (<j>). 
The results of harmonic motions and harmonic ramp motions are plotted in Figs. 3 for 

the delta wing and Figs. 4 for the F-18 model respectively. In the harmonic ramp 
motions, all responses eventually approach the static values corresponding to the angle 
of attack when the motion stops. 

(2)Constant rate pitching motion: 

This is used to compare with test data at the same pitch rate. 

Two special treatments about this type of motion should be mentioned. First of all, an 
arbitrary motion has to be represented locally by a cosine function in order to utilize the 

results of harmonic modeling. At a certain time of arbitrary motion, a and a can be 

described by the cosine and sine functions as 

ay - a /77 +a 0 cos(kt'+<j>) » a m +a (3a) 

a — a 0 ksin(kt' + (i)) (3b) 
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By knowing the harmonic model’s mean angle of attack^) and amplitude(a 0 ), an 

equivalent reduced frequency k and an equivalent phase angle <f> at that given time t’ can 
be solved by Newton’s method. Near the two ends of a harmonic model’s a range, for 
example 2.5 deg. or 62.5 deg. in the present test model, the equivalent k tends to be high 

because of large a and/or a. From the experience with constant rate pitching motions, 

it was found that an unreasonably extrapolated high value of k at a starting point would 
lead to an unacceptable result in simulation. So, one of the variables, a m and a 0 , must 
be treated as an unknown, instead of k, when the extrapolated k-value is greater than a 
given allowable value k max . Through a series of tests, the mean angle of attack a m was 
chosen as the other unknown in case the extrapolated k-value exceeds a given allowable 

Is 

^max* 

Secondly, as the constant-rate pitch-down motions start at a high a, the time integration 
should start from a static value by setting Y -> °° to the first term of Eq. (2). 

The results of the constant rate pitching motions are presented in Figs.5 for the delta 

wing and Figs.6 for the F-18 model respectively. As expected, all results for pitch-up 
motions are well predicted except in the region near the starting point, in particular for C D . 
The reason for this is again the equivalent frequency being too high. Investigation is 
underway to find an improved solution. On the other hand, some of the results for pitch- 
down motions are not as good as those for pitch-up motions. The discrepancies in 
predicting pitch-down motions can be attributed to the poor numerical modeling of 
harmonic motions on down-strokes, especially at low and moderate reduced frequencies. 
The latter is, in turn, caused by the non-smooth variation of the response from one 
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frequency to the other, as can be seen from Figure 7. 

3. Work Underway 

Additional investigation is underway to compare the computed response with data 
when the mean angle of attack and amplitude are significantly different from those used 
in building up the model. 
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Figure 1 Analysis of 70-deg. Delta Wing Dynamic Stall Data 
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(a)Lift Data 
Figure 2 Continued 
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(c) Pitching Moment Data 
Figure 2 Continued 
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(c)Pitching Moment Data 
Figure 2 Concluded 
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Figure 3 Delta Wing Responses by Time Integration for Harmonic Motion 

and Harmonic Ramp Motion 
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(a)Lift Data 

Figure 4 F-18 Responses by Time Integration for Harmonic Motion 
and Harmonic Ramp Motion 
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Figure 5 Delta Wing Responses by Time Integration 
for Constant Rate Pitching Motion 
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Figure 5 Continued 





Figure 5 Continued 




0 10 20 30 40 50 60 70 

a 

(b)Drag Data 
Figure 5 Continued 




-0.5 

0 10 20 


30 40 50 60 70 

a 

(b)Drag Data 


Figure 5 Continued 









































(a)Lift Data 

Figure 6 F-18 Responses by Time Integration for 
Constant Rate Pitching Motion 
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(a) Lift Data 
Figure 6 Continued 
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Fourier Functional Analysis 
for 


Unsteady Aerodynamic Modeling 
Suei Chin and C. Edward Lan 
The University of Kansas 
Lawrence, Kansas 66045 
Abstract 

A method based on Fourier analysis is developed to analyze the force and moment 
data obtained in large-amplitude forced oscillation tests at high angles of attack. The 
aerodynamic models for lift, drag and pitching moment coefficients are built up from a set 
of aerodynamic responses to harmonic motions at different frequencies. The final 
expressions for the models involve time integrals of the indicial type. Results from linear 
two- and three-dimensional unsteady aerodynamic theories as well as test data for a 70- 
deg delta wing are used to verify the models. It is shown that the present modeling 
method is accurate in producing the aerodynamic responses to harmonic motions and the 
ramp-type motions. 

Nomenclature 

coefficient of cosine Fourier series 
coefficient of sine Fourier series 

average value of constant terms in the harmonic oscillation responses 
drag coefficient 
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2- D lift coefficient. 

3- D lift coefficient. 

variation of lift coefficient with respect to angle of attack 

variation of lift coefficient with respect to pitch rate 

pitching moment coefficient 

constants associated with the zero-lag response 

constants in amplitude functions 

imaginary part of a complex number 

index 

reduced frequency (=oofi/v 00 ) 

Mach number 

index for reduced frequency. Also index for the coefficients in Pad6 approximants 

number of frequencies 

Pad6 approximants 

coefficients in Pade approximants 

pitch rate in rad/sec 

time 

nondimensional time (=^^4) 

free stream velocity 

variation in angle of attack (=a Q coskt’) 

= v a 

amplitude of angle-of- attack variation 
mean angle of attack 


& time rate of change in angle of attack 


2 


{ reference length 

x dummy time integration variable 

£ running variable in time 

0 = kt’ 

Introduction 

Due to the requirement of increased performance and maneuverability, the flight 
envelope of a modem fighter is frequently extended to the high angle-of-attack regime. 
Vehicles maneuvering in this regime are subjected to nonlinear aerodynamic loads. The 
nonlinearities are due mainly to three-dimensional separated flow and concentrated 
vortex flow that occur at large angles of attack. Accurate prediction of these nonlinear 
airloads is of great importance in the analysis of a vehicle’s flight motion and in the 
design of its flight control system. As Tobak and Schiff mentioned in ref. 1, the main 
difficulty in determining the relationship between the instantaneous aerodynamic load on 
a maneuvering vehicle and the motion variables is that this relationship is determined 
not only by the instantaneous values of motion variables but also by all of the prior states 
of the motion up to the current state. With the advanced computing techniques, one 
straightforward way to solve this problem is to solve the flow-field and the flight dynamic 
equations simultaneously. However, this is obviously a very costly approach. In 
particular, at high angles of attack the aerodynamic loads depend nonlinearly on the 
motion variables. Under such conditions, even if the vehicles start from closely similar 
initial conditions, they may experience widely varying motion histories. Thus, a 
satisfactory evaluation of the performance envelope of an aircraft may require a large 
number of coupled computations, one for each change in initial conditions. Furthermore, 
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since the motion and the aerodynamic response are linked together in this approach, 
there can be no reutilization of the previously obtained aerodynamic reactions. 

To avoid the disadvantage of solving the coupled flow- field equations and aircraft’s 
motion equations, an alternate approach is to use a mathematical model to describe the 
steady and unsteady aerodynamics for the aircraft’s equations of motion. Ideally, with a 
mathematical model, an evaluation of the aerodynamic terms specified by the model 
would be required only once. The specified model can be reutilized to solve the aircraft’s 
equations of motion over a range of motion variables and flight conditions. 

In the classical linear potential flow theory 2,3 , researchers in the field of 
aeroelasticity used the Fourier transform to relate the aerodynamic response of step 
change in angle of attack of a wing to that of harmonic oscillatory motions. The transient 
aerodynamic reaction to a step change is called the "indicial function" and has been 
calculated for several classes of isolated wings 2 ' 5 . By a suitable superposition 6 of these 
results, the aerodynamic forces and moments induced in any maneuvers can be studied 2,3 
. Tobak applied the indicial function concept to analyze the motions of wings and wing- 
tail combinations 7 . Later, based on a consideration of functional, Tobak and his 
coworkers 1,8 extended the concept of indicial function into the nonlinear aerodynamic 
regimes, even with aerodynamic bifurcations 9 . The simplest nonlinear aerodynamic 
model proposed in ref. 1 has been applied by several authors 10 ' 14 to perform the analysis. 
However, that simplest model is accurate only to the first order of frequency. It needs to 
be improved for a more general response. 

Aerodynamic forces and moments acting on a rapidly maneuvering aircraft are, in 
general, nonlinear functions of motion variables, their time rate of change, and the history 
of maneuvering. How these unsteady aerodynamic forces and moments may be 
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represented in a form suitable for flight dynamic simulation becomes uncertain, in 
particular at high angles of attack. For a certain type of nonlinearities produced in an 
experiment with small-amplitude oscillation, the analysis has been accomplished by 
separating the time-history data into in-phase and out- phase components 15 . When large- 
amplitude forced oscillations are employed in the wind-tunnel testing at a large mean 
angle of attack, the aerodynamic phenomena may involve dynamic stall and/or strong 
vortex flow, with or without vortex breakdown. In this case, higher harmonic components 
in the aerodynamic response are expected to exist 16 and the phenomenon of aerodynamic 
lag would be important. Therefore, a more general modeling technique is needed. 

In this paper, a numerical method will be developed to analyze the nonlinear and 
time-dependent aerodynamic response to establish the generalized indicial function in 
terms of motion variables and their time rates of change. 

Theoretical Development 

In existing flight simulation, two common ways are used to treat high-angle-of 
attack aerodynamics. One is to use tabulated quasi-steady data and the other is to use 
a local linearized model which form a piecewise continuous fit of the nonlinear response 16 . 
In the present approach, a formula involving time integration will be developed. 

Based on functional analysis, Tobak and Schiff 1 developed a fundamental 
formulation of aerodynamic response for arbitrary motion. Summing incremental 
responses to small step changes of a and qfi/v^ at time x, an integral form for C L at time 
t is obtained 


C L (t) - C L (0) + J^ClJcc^), q ($); t, x] d x 

+ — C fc C L [«(£), q ($); t, t] 4s. d x 

Voo q dt 


( 1 ) 
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where £ is a running variable in time over the interval 0 to x, { is a reference length and 
V M is the freestream velocity. 

To have practical applications, the functional integral form needs to be simplified. 

By assuming that a and q are analytical functions in the neighborhood of E=x, variables 

a and q can be expanded by their Taylor series at £=x. The indicial responses C T , for 

a 

example can be expressed as 

ClJoi^), q(^); t, x] - CLjt, x; a(x), &(x), — , 

q(x), 4(x), ] ^ 

If only the first two coefficients are retained in the above Taylor series expansion, the 

integral form of eq. (1) becomes 


C L (t)-C L (0) + J^ t C La [t,x;a(x),<i(x),q(x),<i(x)]^dx 

+ — C t CL n [t,x ;a(x),&(x),q(x),q(x)]^idx 
v M -^ q dx 


(3) 


Eq. (3) is applicable to the study of rapidly varying maneuvers, where hysteresis 
phenomena are known to exist. However, it is difficult to implement eq. (3). By assuming 
a slowly varying motion, Tobak and Schiff neglected the dependence of the indicial 
response on a and q . By further assuming that the indicial response is a function of the 
elapsed time t-x instead of t and x separately, a simplified expression of eq. (1) can be 
written as 


C L (t)-C L (0) + J^ t C Lo [t-T;o(T),q(T)]i^2dT 

Although the form of eq. (4) represents a great simplification over that of eq. (1), the 
equation still includes the full linear form as a special case. 

Jenkins 19 applied a local Taylor expansion to indicial response Cl and used that 
Taylor expansion form to fit numerical indicial responses calculated from a program called 
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NLWAKE. By substituting Ct into eq. (4), Jenkins was able to predict the oscillating 

a 

motion for airfoil at low frequencies. 

In the present investigation, the hysteresis effect is included and the assumption of 
low frequencies will be removed. Therefore, a form between eq. (3) and eq. (4) is written 
as 

ClW-ClCOI+DCt [t-t; oc(t), &(t), q(x), q(r)] ^Idx 

1/0 dT (5) 

+— [t-t; a(t), <i(t), q(t), q(t)] i^lldt 

v M ^ q dt 

In wind-tunnel testing, the q effect cannot be separated from that of a . Since the 
method developed in this study will be used first to analyze the wind tunnel data, a will 
be used instead of q in the following investigation . The effect of bt (i.e. q ) is included in 
the response without aerodynamic lag, such as the virtual mass effect in incompressible 
flow. Since the zero-lag response does not involve the aerodynamic lag, it is removed out 
of the time integral. Then eq. (5) is rewritten for the present study as 
CL(t)-CL(0)+ zero-lag response 

<*m. «*)] ^1* (6) 

o(x), «t)] 

v^ a dt 

The main objective in the present investigation is to find a suitable form for the integrand 
of eq. (6). The basic building blocks of the present method are a set of aerodynamic 
responses to harmonic motions at different frequencies. These responses serve as a 
linearly independent set of functions upon which the response to an arbitrary motion can 
be built. 

In the linear theory 2 ’**, the aerodynamic response can be separated into a product 
of an amplitude function and a phase function in harmonic motion. The amplitude 
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function depends on motion variables and their time rate of change. On the other hand, 
the phase function is a function of frequency and accounts for any phase lag between the 
response and the excitation. In a two-dimensional linear theory, the phase function is 
given by Theodorsen’s circulation function 2 ’^. After response is obtained at different 
frequencies with the same amplitude in harmonic oscillation, the phase function can be 
determined numerically. After use of reciprocal relations 20 , the indicial function can be 
defined by numerical means. This approach has been used for numerical determination of 
indicial lift for plunging airfoils ^ and for plunging wings 21 . 

The method for the linear theory is generalized as follows. Instead of assuming 
that the aerodynamic response is a product of an amplitude function and a phase 
function, it is taken to be a sum of the products of amplitude functions and phase 
functions in harmonic motion; i.e., 

CL-Co+Hamplitude function) j * (phase function) j (7) 

j 

In the linear theory, j equals 1 in the equation. To determine the form of the 
amplitude functions as functions of a(t) and a (t), and the phase functions, a functional 
analysis is needed. A practical method for this purpose is the Fourier analysis of forced 
oscillation data. The motion is assumed to be of the form: 


ai-ajn+aocosfkt^ 

(8a) 

a-aoCosCktO 

(8b) 

A-C-aoklsinCkt^ 

(8c) 


where k is the reduced frequency, t’ is the nondimensionalized time, is the mean angle 
of attack and a Q is the amplitude of angle-of-attack variation. The first step is to Fourier- 
analyze the response over one period. Let 
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C L- Aq +A]COS 0 + A2Cos20 + A3COS30 

+B]sin0+B2sin20+B3sin30 (9) 

+ . • • 

91 99 

From the past experience ’ , it was found that Pad6 approximants provide an accurate 
approximation of the theoretical phase function. Therefore, Padd approximants will be 
used in the present model as phase functions. Following the classical airfoil theory, the 
analysis is best performed in complex algebra. For this purpose, eq. (9) (or the 
experimental oscillatory results) is rewritten in a complex form, as follows: 


CL-Ao^Ai-iB^ei^+CAa-iBale 12 ^ (1( 

+ (A 3 -iB 3 )e i3kt/ + . . . 

It should be kept in mind that only the real part of the response has a physical meaning. 
The reason to use the complex form is to benefit from the mathematical convenience of 
the e* kt notation. If a is rewritten as 

ilrf’ 

a = a Q e 1KX 

and 

a = (ia Q k) e lkt 

then the classical airfoil theory suggests that the response can be put in the following 
form involving the products of amplitude functions and phase functions as 


Cl s Ao(k) 

+ + E21& + Ci * (Hna + H21& ) * (1 - PDi) 

+ Ei2*2 + ^22^2 + C2 * (Hi2<x 2 + H22CX& + H32& 2 ) 

* (1 - pd 2 ) 

+ E 13 & 3 + E23&3 + C 3 * (H 13 a 3 + H 2 3a 2 & + H 33 (x<i 2 + H43& 3 ) 

* (1 - PD 3 ) 


( 11 ) 
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where PD’s are Padd approximants of order 2 and are defined as 

p _ (ik) 2 : p 2j gk) (12) 

P 3j (ik) 2 + (ik) + P 4j 

ctj + E 21 & j etc. are the zero-lag response. The variables dj and are defined as 
d J = ika 0 J e y kt ’ 

and 

a ■ = -k 2 e 1 * 1 ’ 

to be consistent with higher order terms. When j=l in the above equations, a 1 = a and 
a 1 = a . In addition, H 21 , H 22 , H 23 , etc., are related to the pitch-rate effect. It should be 
noted that those terms inside the parentheses following Cp C 2 , C 3 , such as (H-^a + 

H 21 a ), represent the magnitude (or amplitude) and (1 - PDj) represents the unsteady 
aerodynamic lag (or phase) in response. Therefore, the present assumed form for 
aerodynamic modeling encompasses the classical linear theory and is capable of 
representing a complete set of harmonic-oscillatory data with different frequencies in one 
expression. It should be noted that in eq. (11), the contribution to each mode is summed 
in complex form. The response in time domain is given by the real part, similar to 
obtaining eq. (9) from eq. (10). 

Cj are the reference values used to normalize the lift given by Aj - i Bj in the least 
squared-error method, j is the index consistent with the exponent of the exponential 
terms in eq. (10). For example if the j’s term in eq. (11) represents the coefficient of e , 
then j is 1. If the j’s term in eq. (11) represents the coefficient of e i2kt then j is 2, etc. 

The first term, A^k), in eq. (11) is a constant term, supposedly a function of frequency. 
From available experimental data for a delta wing 22 , Ag(k) can be assumed to be constant 
approximately. The unknown coefficients P^, P 2 j, P 3j - and P 4 j are calculated from the 
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least squared-error method. E^, E 21 , H-^, H 12 , etc., are obtained separately by 
minimizing the sum of squares of errors. This is equivalent to a two-level optimization 
method to determine the unknowns in eq. (11). That is, E, H, etc., are assumed first. 

Then Pjj, etc., are determined by minimizing the sum of squared errors. The values of 
En, ^11* e ^ C-> 33:8 var ^ e< ^ nex t so that the sum of squared errors is minimized based on a 
gradient method. It was found that this approach is more effective in determining a 
global minimum solution for the unknowns than a straightforward optimization (one 
level) method because of nonlinearity in the unknowns in this optimization problem. It 
should be noted that in the literature the phase function has been typically determined by 
the response to plunging motions, not pitching motions. Therefore, those terms associated 


with a in eq. (11) do not appear. This would very much simplify the mathematics of 
determining the Pad6 approximants. The details of the present method are discussed in 


the following. 
Least-Square Method 


By choosing proper values of E^, Hu, H 12 , etc., in eq. (11), the corresponding Aj - 
i Bj term in eq. (10) is then divided by the amplitude function. The result will appear as 


y. + iw- = 1 - A J - iB j - E lj & ~ E 2j (-k 2 ) 

J J (amplitude function); 

(I 3 ) 

Plj (ik) 2 + P2j (ik) 

P 3j (ik) 2 7 (ik) 7 P 4j 

If both sides of eq. (13) are multiplied by the denominator of the Padd approximant and 


separated into real and imaginary parts, then 
Re * V 2 - P 3 j v j k2 + P 4 j Vj - W) k = 0 

and 

Im.P 2j k + P 3j W j k 2 .P 4j W r V j k = 0 


(14a) 

(14b) 
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The sum of squared errors is defined as 

Err s E Re(ki) 2 + E Im(kj) 2 (15) 

By equating the first derivatives of squared errors (eq. 15) with respect to variables 

Pjj, P 2 j , P 3 j and P 4j - to zero, the unknown coefficients P^, P 2 j, P 3 j and P 4j - can be 

determined from 

Ekf 0 -EVjkf EVjkf 

0 Ekf EWjkf -EWjkj 

-EVjkf EWjkf E(vfkf+wfkf) -E(vfkf+wfkf 

EVjkf -EWjkj E(V?kf+wfkf) E(V?+wf) 

where i varies over the range of input frequencies, and the mode subscript j on V and W 
has been omitted. 

Gradient Method 

After the unknown coefficients P-y, P 2 j, P 3j - and P 4 j have been found, a one- 
dimensional gradient method is used to find E and H values which will make the sum of 
the squared errors minimum. The E or H value is perturbed first by a small amount AE 
or AH to find the gradient of the sum of squared errors. If the gradient tends to reduce 
the error, then the E or H value is perturbed further until several iterations has been 
reached (it is set to be 5 iterations in the current program). After that, the same 
procedure is applied to other E or H. Then the whole procedure is repeated again for 
several iterations. 

Indicial Formulation 

To express the aerodynamic response in time domain (eq. 6), the phase function, as 
represented by the Padd approximants, is inverted from frequency to time domains by 
inverse Fourier transform. The Pad6 approximants are first factored as follows: 
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( 17 ) 


Plj (ik) 2 + P2j (ik) 1 _ ik axj ik a2j 
P 3j (ik) 2 + ik + P 4j ik + a 3j ik + a 4j 

and then it is inverted based on a step input 2 to be 


f-[l - Pl € k _— M— ] e «'d (jk) 

2tc xk P 3j (ik) 2 + ik + P 4j 

- 1 - ay e” a3 i fc - a2j e" a4 i fc 

The final form for the aerodynamic response in time domain for arbitrary motion is 
therefore given by 


m 


C L( fc/ ) - CLindiciaitt'-T, + C ave + E (Eij&j + E 2 jixj) 

j-1 

. E ft' d(a - f) i.(l-a li e- a ^ t, - l) -a 2j e-^ t '- 1> ) *£> dx 
j- 1 ^ 0 da J J dx 

+ J.E ft' li!£l.(l-a u e- aa i (t '- t) -a 2i e- a #'- ,) ) i dx 
Veoj-l*' 0 d 6 J J dx 


(19) 


where the amplitude functions (a.f.) are given by those H-terms in eq. (11) and C ave is the 
average of Aq and is a function of a^. It should be noted that a in eq. (19) denotes a 
perturbation from o^. The first term in eq. (19) is the amplitude of C L (eq. 11) when a is 
abruptly changed to a(0) at x = 0 and represents an initial value in the indicial lift 
formulation (see eqs. 5.370 and 5.382 of Ref. 2). Again, each mode is evaluated in 
complex form and the real part of the result is taken as the response in time domain. 

To perform the time integration in eq. (19), the 3-point Simpson rule is used in the 
present method. Since the amplitude functions are determined in the frequency domain 
using complex algebra, for an arbitrary motion an equivalent frequency k and phase angle 
<t> at x must be obtained by matching the given a-^ and a ^ as follows: 
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a l(x )-a m +a o cos(kx +<j> ) 
(3c i — a 0 ksin(kx+<t>) 


(20a) 


(20b) 

Eqs. (20) are solved by Newton’s method. It should be noted that k and <$> are needed 
merely to simulate an equivalent harmonic motion in the present formulation. The 
resulting k and <j> are then used to determine the magnitude of the amplitude function 
using complex algebra. 

Results and Discussions 

Because appropriate high-alpha experimental data to apply the present modeling 
method are limited, the present method will first be tested with linear theoretical results. 
Linear Results 

Several cases in the two-dimensional and three-dimensional linear flow at different 
Mach numbers have been studied to verify the present method of aerodynamic modeling. 
The oscillatory 2-D results are computed from a 2-D unsteady quasi-vortex-lattice method 
(QVLM) 24 as input data for the current model. Through numerical experimentation, it is 
found that six frequencies are needed to have accurate results. In the present model (eq. 
11), only the coefficients E^, E 21 , H^, H 21 and are not zero. Three Mach numbers, 
0, 0.2 and 0.4 are employed. The results at M=0.4 are presented in Fig. 1. It is seen that 
the present modeling method is very accurate for harmonic motion. The modeling method 
is further verified with a 70-degree delta wing which oscillates from 0 to 20 deg. in angle 
of attack about the mid-root chord: 

a 1 = 0.17453 + 0.17453 cos kt’ (in radian) 

This means that the mean angle of attack is 10 deg. (0.17453 radian) and the amplitude 
of the oscillation is also 10 deg. (0.17453 radian). The aerodynamic responses are 

AC 

calculated from a 3-D unsteady QVLM code . Through numerical experimentation, it is 
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found that the responses at low frequencies do not change significantly, so that accuracy 
in modeling would be reduced. To have more accurate approximation, high frequencies’ 
responses are needed. Seven reduced frequencies (k = 0.01, 0.1, 0.2, 0.6, 1.0, 2.0, 2.5) are 
used in this 3-D attached flow case. The results for at M=0.4 from modeling are 
plotted in Fig. 2 and show very good agreement again with results from the 3-D unsteady 
QVLM program. 

Nonlinear Results 

oo 

The forced-oscillation test data for a 70-deg delta wing in pitching oscillation are 
used to validate the present nonlinear aerodynamic model. The angle of attack which 
describes the pitching motion is given as 

<Xj_ = 27.5 - 27.5 cos kt’ (in degree) 

which means the delta wing oscillates from 0 to 55 deg. in angle of attack and then back 
to 0 deg. for one cycle. The pitching center is at 57 % of the root chord. The reduced 
frequency k is nondimensionalized based on wing’s root chord. Five sets of data 
corresponding to 5 different frequencies are available and they are used as the input data 
to calculate the coefficients for the current aerodynamic model with 5 Fourier terms. The 
lift coefficients obtained from the aerodynamic model (eq. 11) are compared with the 
original test data in Fig. 3 with good agreement. Expressions for C D and C m similar to 
eq. (11) are obtained with the same procedures as those used for C L . The modeled 
harmonic results are compared with data in Figs. 4 and 5 for Cj-j and C m , respectively. 
Again, the good agreement indicates that the present aerodynamic model is accurate in 
representing the experimental harmonic data. The coefficients for C L are tabulated in 
Table 1. 

To check the validity of indidal formulation (eq. 19) for the present nonlinear 
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response, two oscillatory cases in Fig. 3 will be used. That is, by assuming oscillatory 
motion in eq. (19), the time-integrated lift response should agree with the forced- 
oscillation results. The lift coefficient by integrating eq. (19) for the same 70-deg. 
oscillating delta wing with k=0.098 is plotted in Fig. 6. Compared with the results from 
aerodynamic modeling, the integrated lift shows good agreement. 

To verify the aerodynamic models further, aerodynamic responses to harmonic 
ramp motions for a 70-deg delta wing reported in ref. 26 will be employed. The ramp 
motions start from a = 0 to 35, 45 and 55 deg. In the present calculation based on eq. 

(19), the same harmonic data shown in Figs. 3 - 5 and reported in ref. 23 are used. That 
is, the data are based on harmonic motions from a = 0 to 55 deg. The results for C L are 
presented in Fig. 7 for a reduced frequency of 0.0714. It is seen that the present 
aerodynamic model is fairly accurate if the harmonic ramp motion is from a = 0 to 55 deg. 
However, the final C L is overpredicted if the ramp motion stops at an a less than 55 deg, 
even though the peak Cl is still well predicted. A possible reason for this is that the 
harmonic data based on a = 0 - 55 deg. contain dynamic effect on vortex-breakdown 
characteristics at a < 55 deg. Therefore, the results for Cl at a final steady a = 35 or 45 
deg. should be higher. 

The corresponding drag and pitching moment coefficients at one reduced frequency 
are presented in Figs. 8 and 9, respectively. The drag coefficient is not as well predicted 
in ramp motions as in harmonic motions (see Fig. 4). It is not known whether this is 
caused by differences in the test models and test Reynolds numbers. The test model for 
the harmonic motions (ref. 23) has two-sided chamfered leading edges with a thickness of 
0.5 inch at a Reynolds number of 1.64xl0 6 based on the root chord. The model for the 
ramp motions (ref. 26) is chamfered only on the lower surface of the leading edge and has 
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a thickness of 0.25 inch, and tested at a Reynolds number of 1.54x10®. The pitching 
moment coefficient appears to be well predicted except at small time. 

To illustrate the present aerodynamic model (eq. 19) for arbitrary motions, a linear 
ramp motion is assumed in the integration. The results are compared with those in a 
harmonic ramp motion in Fig. 10. It is seen that the linear ramp motion tends to produce 

higher Cl beyond the peak value because it has a higher value in tic . 


Although verification of the present model was presented only with one set of cc Q 
and cc^, some preliminary results with different a Q and indicate that eq. (19) could 
still produce good results if a new a-range is within the test range used in setting up the 
model. 


Table 1 Model Coefficients for C L for the 70-deg Delta Wing 
c ave = 0.6451 


j 

C i 

HI 

■Mi 

KTTfli 

mm 

wm 

mm 

KM 

KM 

m 

1 

1.0 

-0.389 

1.062 

0.700 






2 

1.0 

0.212 

0.250 

-0.700 

0.500 

0.600 




3 

1.0 

-0.368 

0.118 

-0.970 

0.534 

0.995 

-1.019 



4 

5.0 

0.025 

-0.063 

-0.100 

0.400 

0.400 

1.000 

0.000 


5 

30.0 

0.186 

0.039 

0.096 

0.588 

-0.015 

-0.020 

.0009 

.0007 

j 

wm 

wm 

wm 

wm 

IMj 

■N 

KM 



1 

Egg 

-0.453 

iSRWi 


■BPS] 

-0.646 

-0.037 



2 

4.947 

-1.387 

15.243 

0.001 

-1.437 

1.761 

-0.001 

-0.065 


3 

3.561 

0.653 

4.383 

0.041 

0.866 

-0.054 

-0.053 

-0.175 


4 

24.424 

3.412 

21.343 

0.001 

3.541 

-2.396 

-0.001 

-0.046 


5 

6.127 

1.604 

1.244 

0.025 

1.545 

3.379 

-0.026 

-0.778 
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Conclusions 


A Fourier analysis method was developed to analyze harmonic forced-oscillation 
data at high angles of attack as functions of the angle of attack and its time rate of 
change. The resulting aerodynamic responses at different frequencies are used to build 
up the aerodynamic models involving time integrals of the indicial type. An efficient 
numerical method was also developed to evaluate these time integrals for arbitrary 
motions based on a concept of equivalent harmonic motion. The method was verified by 
first using results from two-dimensional and three-dimensional linear theories. The 
developed models for C^, Cq and C m based on high-alpha data for a 70-deg delta wing in 
harmonic motions showed accurate results in reproducing hysteresis. The aerodynamic 
models are further verified by comparing with test data using ramp-type motions. 
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Figure 1 Unsteady Lift Coefficient for a 2-D Flat Plate Pitching about Midchord at M = 

0.4 

Figure 2 Unsteady Lift Coefficient for a 70-deg Delta Wing Pitching about Mid-root Chord 
at M = 0.4 

Figure 3 Unsteady Lift Coefficient for a 70-deg Delta Wing Pitching about 57% of Root 
Chord at Low Speed and Various Frequencies. Re = 1.64x10® 

Figure 4 Unsteady Drag Coefficient for a 70-deg Delta Wing Pitching about 57% of Root 
Chord at Low Speed and Various Frequencies. Re = 1.64x10® 

Figure 5 Unsteady Pitching Moment Coefficient for a 70-deg Delta Wing Pitching about 

57% of Root Chord at Low Speed and Various Frequencies. Moment Center at 25% 
of Root Chord. Re = 1.64x10® 

Figure 6 Unsteady Lift Coefficient from Numerical Modeling and Indicial Time 

Integration for a 70-deg Delta Wing in Harmonic Pitching Oscillation about 57% of 
Root Chord at Low Speed. Re = 1.64x10® and k = 0.098 

Figure 7 Unsteady Lift Coefficient from Indicial Lift Model and Experiment for a 70-deg 
Delta Wing in Harmonic Ramp Motion at Low Speed. Re = 1.54x10® and k = 
0.0714 

Figure 8 Unsteady Drag Coefficient from Indicial Drag Model and Experiment for a 70- 
deg Delta Wing in Harmonic Ramp Motion at Low Speed. Re = 1.54x10® and k = 
0.0714 

Figure 9 Unsteady Pitching Moment Coefficient from Indicial Pitching Model and 

Experiment for a 70-deg Delta Wing in Harmonic Ramp Motion at Low Speed. 
Moment Center at 25% of Root Chord. Re = 1.54x10® and k = 0.0714 

Figure 10 Unsteady Lift Coefficient from Indicial Lift Model in Harmonic and Linear 
Ramp Motions for a 70-deg Delta Wing at Low Speed. Re = 1.54x10® and k = 0.0926 
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